SEA rsfcMRI - Longitudinal Descriptives

if(fslong){
  cat('\n\n# FreeSurfer Longitudinal Registration\n\n')
} else {
  cat('\n\n# Subject-session Registration\n\n')
}

FreeSurfer Longitudinal Registration

Load the data

Starting with the Schaefer 400-parcel, 7 network atlas.

Retain just within-network connectivity. Also, Fisher-z transform the correlations for analyses. One small detail that is important here is that we keep network connectivity for same-hemisphere parcels only.

library(data.table)
library(dplyr)
library(tidyr)

#set variables for `collect_data`
if(fslong){
  fname <- 'sea_rsfc_fslong_schaefer400x7.RDS'
  data_dir <- '/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/SEA_REST/derivatives/fmriprep-20.1.1/xcpengine-default/'
  subs <- sprintf('sub-10%02d', c(1:16,18:31))
  sess <- sprintf('%02d', 1:10)
  subses <- expand.grid(sess, subs)
  files <- paste0(data_dir, '/', subses[,2], subses[,1], '/fcon/schaefer400x7/', subses[,2], subses[,1], '_schaefer400x7.net')
  file_id <- paste0(subses[,2], '_', subses[,1])  
} else {
  fname <- 'sea_rsfc_schaefer400x7.RDS'
  data_dir <- '/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/SEA_BIDS/derivatives/1.4.1-final/xcpengine-default'
  subs <- sprintf('sub-10%02d', c(1:16,18:31))
  sess <- sprintf('ses-%02d', 1:10)
  subses <- expand.grid(sess, subs)
  files <- paste0(data_dir, '/', subses[,2], '/', subses[,1], '/fcon/schaefer400x7/', subses[,2], '_', subses[,1], '_schaefer400x7.net')
  file_id <- paste0(subses[,2], '_', subses[,1])
}

library(data.table)
library(tidyr)
if(is.na(ncores <- as.numeric(Sys.getenv('SLURM_CPUS_ON_NODE')))){
  ncores <- 3
  NOTSLURM <- TRUE
  message('No environment variable specifying number of cores. Setting to ', ncores, '.')
} else {
  NOTSLURM <- FALSE
  message('Number of cores set to ', ncores)
}
setDTthreads(ncores)

if(file.exists(fname)){
  adt_labels <- readRDS(fname)
  schaefer_400_7_net_ids <- readRDS('schaefer_ids.RDS')
} else {
  message('Creating file list...')
  files_list <- lapply(files, function(x) ifelse(file.exists(x), x, ''))
  
  message('Reading rsfc files...')
  if(ncores > 1){
    message('Using ', ncores, ' cores to read files.')
    cl <- parallel::makePSOCKcluster(max(c(ncores - 1), 1))
    
    adt_list <- unlist(parallel::parLapply(cl = cl, split(files, 1:ncores), function(part){
      lapply(part, function(f){
        if(file.exists(f)){
          data.table::fread(f, skip = 2, col.names = c('node1', 'node2', 'r'))
        } else {
          data.table::data.table()
        }
      })
    }), recursive = FALSE)
    try(parallel::stopCluster(cl))
  } else {
    adt_list <- lapply(files, function(f){
      if(file.exists(f)){
        data.table::fread(f, skip = 2, col.names = c('node1', 'node2', 'r'))
      } else {
        data.table::data.table()
      }
    }) 
  }
  names(adt_list) <- file_id
  
  message('Combining data, assigning labels...')
  adt <- rbindlist(adt_list, idcol = 'id')
  
  #https://github.com/ThomasYeoLab/CBIG/tree/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/Parcellations/MNI
  schaefer_400_7_net_ids <- fread('Schaefer2018_400Parcels_7Networks_order.txt', 
                                  col.names = c('node1', 'network1', 'V1', 'V2', 'V3', 'V4'))
  schaefer_400_7_net_ids[, c('V1', 'V2', 'V3', 'V4') := NULL]
  schaefer_400_7_net_ids <- schaefer_400_7_net_ids %>% 
    tidyr::extract(network1, c('network1', 'roi1'), '7Networks_([RL]H_.*)_(\\d+)')
  
  setkey(adt, node1)
  setkey(schaefer_400_7_net_ids, node1)
  adt_labels1 <- schaefer_400_7_net_ids[adt]
  setnames(schaefer_400_7_net_ids, c('node2', 'network2', 'roi2'))
  setkey(adt_labels1, node2)
  setkey(schaefer_400_7_net_ids, node2)
  adt_labels <- schaefer_400_7_net_ids[adt_labels1]
  adt_labels[, c('id', 'sess') := tstrsplit(id, '_', fixed = TRUE)]
  message('Saving RDS file to: ', fname)
  saveRDS(adt_labels, fname)
  message('Saving Schaefer ID data table to schaefer_ids.rds')
  saveRDS(schaefer_400_7_net_ids, 'schaefer_ids.RDS')
}

sea_schaefer400_7_withinnet <- copy(adt_labels)

sub_regex <- '[LR]H_[A-Za-z]+_*(.*)'
net_regex <- '[LR]H_([A-Za-z]+)_*(?:.*)'
hem_regex <- '([LR]H)_[A-Za-z]+_*(?:.*)'
sea_schaefer400_7_withinnet_nets <- distinct(sea_schaefer400_7_withinnet, network1)
sea_schaefer400_7_withinnet_nets[, sub1 := gsub(sub_regex, '\\1', network1)]
sea_schaefer400_7_withinnet_nets[, net1 := gsub(net_regex, '\\1', network1)]
sea_schaefer400_7_withinnet_nets[, hem1 := gsub(hem_regex, '\\1', network1)]
setkey(sea_schaefer400_7_withinnet, network1)
setkey(sea_schaefer400_7_withinnet_nets, network1)
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet_nets[sea_schaefer400_7_withinnet]

setnames(sea_schaefer400_7_withinnet_nets, c('network2', 'sub2', 'net2', 'hem2'))
setkey(sea_schaefer400_7_withinnet, network2)
setkey(sea_schaefer400_7_withinnet_nets, network2)
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet_nets[sea_schaefer400_7_withinnet]

sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet[net1 == net2 & hem1 == hem2]
sea_schaefer400_7_withinnet[, c('network1', 'network2') := NULL]
sea_schaefer400_7_withinnet[, z := atanh(r)]
sea_schaefer400_7_withinnet[, H := dplyr::case_when(hem1 == 'RH' ~ -1,
                                                    hem1 == 'LH' ~ 1)]
if(!dim(distinct(sea_schaefer400_7_withinnet, net1, net2, hem1, hem2))[[1]] == 14){
  stop('Incorrect number of networks')
}

if(FALSE){
  adt_labels_old <- readRDS('sea_rsfc_schaefer400x7.RDS')
  missing_files_csv <- data.table(file = files, does_not_exist = files_list == '')
  View(missing_files_csv[does_not_exist == TRUE])
  readr::write_csv(data.table(file = files, does_not_exist = files_list == ''), 
                   '~/sea_rsfc-missing_fcon_output.csv')
  
  a <- data.table(file = files, does_not_exist = files_list == '')
  a[, fcon_missing := lapply(gsub('(.*/fcon)/.*', '\\1', file), function(x) !file.exists(x))]
  a[, missmatch := does_not_exist != fcon_missing ]
  a[(missmatch), file]
}

Check No. of RS Features

Using the number nodes are in each network allows computation of the expected number of rsfc features, which allows a comparison to the observed number of features to find instances where certain connections are missing.

#Check to make sure we're getting the number of connectivity features per network we expect
setkey(sea_schaefer400_7_withinnet_nets, 'network2')
setkey(schaefer_400_7_net_ids, 'network2')

nodes_per_net <- sea_schaefer400_7_withinnet_nets[schaefer_400_7_net_ids][, .N, by = c('net2', 'hem2')]
nodes_per_net[, expected_Nfeatures := N * (N - 1) / 2]

connectivity_feature_check <- sea_schaefer400_7_withinnet[, .(Nfeatures = .N), by = c('net2', 'hem2', 'id', 'sess')][nodes_per_net[, !'N'], on = c('net2', 'hem2')]

setnames(connectivity_feature_check, c('net2', 'hem2'), c('nework', 'hemisphere'))

dt_options <- list(rownames = FALSE,
                   filter = 'top',
                   class = 'cell-border stripe',
                   extensions = 'Buttons', 
                   options = list(dom = 'Bfrtip', buttons = c('csv')))
do.call(DT::datatable, 
        c(list(connectivity_feature_check[Nfeatures != expected_Nfeatures],
               caption = 'Mismatch of observed and expected number of rsfc features'), 
          dt_options))

do.call(DT::datatable, 
        c(list(connectivity_feature_check[Nfeatures != expected_Nfeatures][, .N, by = c('id', 'sess')],
               caption = 'Number of mismatches for subject sessions'),
          dt_options))

Node coverage

Each node in a network should be connected to N - 1 other nodes in the network, where N is the number of total nodes in the network.

Mean connections-per-node

node_check <- sea_schaefer400_7_withinnet[, c('node1', 'node2', 'id', 'sess', 'net2', 'hem2')]
node_check[, idx := 1:.N, by = c('id', 'sess', 'net2', 'hem2')]
node_check[, c('maxnode', 'minnode') := .(max(node1, node2), min(node1, node2)), by = c('node1', 'node2')]
node_check_l <- melt(node_check, measure.vars = c('maxnode', 'minnode'), value.name = 'node')
node_info <- node_check_l[, .N, 
                          by = c('id', 'sess', 'net2', 'node', 'hem2')][, list(total = max(N)), 
                                                                        by = c('node', 'net2', 'hem2')]
net_info <- node_info[, .(total = unique(total)), by = c('net2', 'hem2')]
node_check_subj <- node_check_l[, .N, by = c('sess', 'id', 'node', 'net2', 'hem2')]

node_check_subj <- net_info[node_check_subj, on = c('net2', 'hem2')]
node_check_subj_summary <- node_check_subj[, .(observed = mean(N),
                                               expected = mean(total),
                                               proportion = sum(N) / sum(total)),
                                           by = c('node', 'net2', 'hem2')]
node_check_subj_summary[, network := paste(net2, hem2, sep = '_')]

node_check_subj_summary_l <- melt(node_check_subj_summary, measure.vars = c('observed', 'expected'))
node_check_subj_summary_l[, variable := factor(variable, levels = c('expected', 'observed'))]
ggplot(node_check_subj_summary_l,
       aes(x = node, y = value)) + 
  geom_line(aes(fill = variable, color = variable), position = position_identity()) + 
  scale_color_manual(values = apal[c(1,2)], aesthetics = c('color', 'fill')) + 
  facet_wrap(~ network, ncol = 2, scales = 'free') + 
  jftheme + 
  labs(y = 'Mean number of connections from node', x = 'Node')

Proportion connections-per-node

ggplot(node_check_subj_summary,
       aes(x = node, y = proportion)) + 
  geom_line(color = apal[[4]]) + 
  facet_wrap(~ network, ncol = 2, scales = 'free') + 
  jftheme + 
  labs(y = 'Proportion of expected connections from node', x = 'Node')

Mean across sessions by subject

node_check_sess_summary <- node_check_subj[, .(mean = mean(N),
                                               proportion = sum(N) / sum(total)),
                                           by = c('id', 'net2', 'hem2', 'sess')]
node_check_sess_summary[, sess_num := as.numeric(gsub('ses-(\\d+)', '\\1', sess))]
ggplot(node_check_sess_summary,
       aes(x = sess_num, y = mean)) + 
  geom_line(aes(group = id), color = apal[[1]], alpha = .5) + 
  facet_grid(net2 ~ hem2, scales = 'free') + 
  jftheme + 
  labs(y = 'Mean variability per node', x = 'Session')

Proportion across sessions by subject

node_check_sess_summary <- node_check_subj[, .(mean = mean(N),
                                               proportion = sum(N) / sum(total)),
                                           by = c('id', 'net2', 'hem2', 'sess')]
node_check_sess_summary[, sess_num := as.numeric(gsub('ses-(\\d+)', '\\1', sess))]
ggplot(node_check_sess_summary,
       aes(x = sess_num, y = proportion)) + 
  geom_line(aes(group = id), color = apal[[1]], alpha = .5) + 
  facet_grid(net2 ~ hem2) + 
  jftheme + 
  labs(y = 'Proportion variability per node', x = 'Session')

Some QC

Density of functional connectivity (pearson correlation) for each participant, for each subject, overlayed on the density for all sessions and participants.

#Generate group-level density
agg_mean <- mean(sea_schaefer400_7_withinnet$z)
agg_sd <- sd(sea_schaefer400_7_withinnet$z)
scalez <- (sea_schaefer400_7_withinnet$z - agg_mean) / agg_sd
density_agg <- density(scalez)
sea_schaefer400_densplot_agg <- data.table(x = density_agg$x,
                                           x_on_z = density_agg$x * agg_sd + agg_mean,
                                           y = density_agg$y)
#Generate per-participant-session densities
sea_schaefer400_MSD <- sea_schaefer400_7_withinnet[, list(mean = mean(z),
                                                          sd = sd(z)), 
                                                   by = c('id', 'sess')]
sea_schaefer400_densplot <- sea_schaefer400_MSD[sea_schaefer400_7_withinnet, 
                                                on = c('id', 'sess')]
sea_schaefer400_densplot[, scalez := (z - mean)/sd]
sea_schaefer400_densplot <- 
  sea_schaefer400_MSD[sea_schaefer400_densplot[, list(x = density(scalez)$x,
                                                      y = density(scalez)$y), 
                                               by = c('id', 'sess')],
                      on = c('id', 'sess')][, x_on_z := sd*x + mean]

u_ids <- unique(sea_schaefer400_densplot$id)
for(an_id in u_ids){
  cat(paste0('\n\n## ', an_id, '\n\n'))
  
  hplot <- ggplot(sea_schaefer400_densplot_agg, aes(x = tanh(x_on_z), y = y)) + 
    geom_ribbon(ymin = 0, aes(ymax = y, fill = 'Group', color  = 'Group'), 
                alpha = .5) + 
    geom_ribbon(data = sea_schaefer400_densplot[id == an_id, ], 
                ymin = 0, aes(ymax = y, fill = 'ID', color  = 'ID'), 
                alpha = .5) + 
    scale_fill_manual(aesthetics = c('color','fill'), breaks = c('Group', 'ID'), labels = c('Group', 'Participant'), values = apal[c(2,4)], name = 'Data') +
    facet_wrap(~sess, ncol = 2) + 
    labs(x = 'correlation', y = 'density') + 
    coord_cartesian(xlim = c(-1, 1), ylim = c(0, .5)) + 
    jftheme
  print(hplot)
}

sub-1001

sub-1002

sub-1003

sub-1004

sub-1005

sub-1006

sub-1007

sub-1008

sub-1009

sub-1010

sub-1011

sub-1012

sub-1013

sub-1014

sub-1015

sub-1016

sub-1018

sub-1019

sub-1020

sub-1021

sub-1022

sub-1023

sub-1024

sub-1025

sub-1026

sub-1027

sub-1028

sub-1029

sub-1030

sub-1031

Modeling

We are interested in stability and variability for each network. For each participant, we have multiple sessions, and within each session, we have multiple parcel-parcel pairs that provide information about the within-network connectivity. We can use the fact that we have multiple indicators of within-network connectivity to estimate the between-person variability, as well as the within-person (session-to-session) variability as distinct from measurement error (we assume that each parcel-parcel functional connectivity estimate in a network is an estimate of that network's connectivity)^*.

* This assumption means that we essentially treat differences in the FC between one pair and another as measurement error.

We can compute an ICC that describes both within-person and between-person variability by using a 3 level model. First, I subset the data for a single network. I then build an intercept-only model, allowing the intercept to vary by participant-ID, and by session within each ID. The intercept is effectively the mean across all intrahemispheric parcel-pairs.

\[ \begin{align} z &= \beta_{0jk} + \epsilon_{ijk} \\ \beta_{0jk} &= \pi_{00k} + \nu_{0jk} \\ \pi_{00k} &= \gamma_{000} + \upsilon_{00k} \end{align} \]

\[ \begin{align} z = \gamma_{000} + \upsilon_{00k} + \nu_{0jk} + \epsilon_{ijk} \end{align} \] Where \(z\) is the measure of functional connectivity, and \(i\) indexes observations within each session, \(j\), for each participant \(k\). This means that we are able to estimate variance in per-person deviations (\(\upsilon_{00k}\)) from the network mean-connectivity (\(\gamma_{000}\)), per-session-deviations (\(\nu_{0jk}\)) from those per-person deviations, and the residual not explained by variability over persons or person-sessions (\(\epsilon_{ijk}\)).

Testing these out

I want to make sure that these models are correct, so I'll compare the variance portions, and total variance, between the two models.

library(lme4)
networks <- unique(sea_schaefer400_7_withinnet$net1)
this_net <- networks[[2]]
this_net_dt <- sea_schaefer400_7_withinnet[net1 == this_net] 
fslong_insert <- ifelse(fslong, '_fslong', '')

model_dir <- 'models'
test_model_list <- list(test_2l = list(fn = file.path(model_dir, paste0('test_noh', fslong_insert, '_2l.RDS')),
                                         form = z ~ 1 + (1 | id)),
                          test_3l = list(fn = file.path(model_dir, paste0('test_noh', fslong_insert, '_3l.RDS')),
                                         form = z ~ 1 + (1 | id/sess)))
cl <- parallel::makePSOCKcluster(max(c(ncores - 1), 1))

test_model_fits <- parallel::parLapply(cl = cl, test_model_list, function(model_list, d){
  library(lme4)
  if(!file.exists(model_list[['fn']])){
    f <- model_list[['form']]
    fit <- lmer(f, data = d)
    saveRDS(fit, model_list[['fn']])
  } else {
    fit <- readRDS(model_list[['fn']])
  }
  return(fit)
}, d = this_net_dt)

try(parallel::stopCluster(cl))

three_level <- test_model_fits$test_3l

summary(three_level)

## Linear mixed model fit by REML ['lmerMod']
## Formula: z ~ 1 + (1 | id/sess)
##    Data: d
## 
## REML criterion at convergence: 181397.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7510 -0.6822 -0.0861  0.5938  6.3337 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sess:id  (Intercept) 0.005088 0.07133 
##  id       (Intercept) 0.000000 0.00000 
##  Residual             0.081313 0.28516 
## Number of obs: 548180, groups:  sess:id, 281; id, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.185122   0.004273   43.32
## convergence code: 0
## boundary (singular) fit: see ?isSingular

lll_varcorr <- VarCorr(three_level)

report_df <- data.frame(
  stat = c('id_var', 'sess_var', 'resid_var', 'total_var'),
  three_level = c(lll_varcorr$id.1['(Intercept)','(Intercept)'],
                  lll_varcorr$sess.id.1['(Intercept)','(Intercept)'],
                  sigma(three_level)^2,
                  sum(unlist(lapply(lll_varcorr, diag))) + sigma(three_level)^2))

report_df$three_level_perc <- 
  sprintf('%.1f%%', report_df$three_level/report_df$three_level[[4]]*100)
report_df$three_level_RE_perc <- 
  c(sprintf('%.1f%%', report_df$three_level[1:2]/sum(report_df$three_level[1:2])*100), '-', '-')
knitr::kable(report_df, digits = 5)

stat	three_level	three_level_perc	three_level_RE_perc
id_var	0.08131	94.1%	48.5%
sess_var	0.08640	100.0%	51.5%
resid_var	0.08131	94.1%	-
total_var	0.08640	100.0%	-

Fit for each network

networks <- unique(sea_schaefer400_7_withinnet$net1)

netlist <- lapply(networks, function(this_net){
  return(list(network = this_net, d = sea_schaefer400_7_withinnet[net1 == this_net]))
})

if(NOTSLURM){
  message('Not running as batch job. Will only attempt to read saved model fits.')
}

cl <- parallel::makePSOCKcluster(max(c(ncores - 1, 1)))
model_fits <- parallel::parLapply(cl = cl, netlist, function(anet, fslong, NOTSLURM){
  this_net <- anet[['network']]
  this_net_dt <- anet[['d']]
  model_dir <- 'models'
  fslong_name <- ''
  if(fslong)
    fslong_name <- 'fslong-'
  model_fit_fn <- file.path(model_dir, paste0('lmer_fit-3l-', fslong_name, this_net, '.RDS'))
  library(lme4)
  if(!file.exists(model_fit_fn) & NOTSLURM){
    stop('Model fit has not been saved. Will not estimate. Stopping.')
  }
  if(!file.exists(model_fit_fn)){
    message('Fitting model for ', this_net)
    fit <- lmer(z ~ 1 + (1 | id/sess), data = this_net_dt)
    saveRDS(fit, model_fit_fn)
  } else {
    fit <- readRDS(model_fit_fn)
  }
  return(fit)
}, fslong = fslong,  NOTSLURM = NOTSLURM)
try(parallel::stopCluster(cl))

proportion_variance_tables <- lapply(model_fits, function(model_fit){
  vc <- VarCorr(model_fit)
  id_v <- vc$id['(Intercept)','(Intercept)']
  idsess_v <- vc$`sess:id`['(Intercept)','(Intercept)']
  s_2 <- sigma(model_fit)^2
  total_RE <- sum(unlist(lapply(vc, diag)))
  other_RE <- total_RE - id_v - idsess_v
  total <- total_RE + s_2
  rez <- data.frame(source = rep(c('ID', 'ID/Session', 'Other RE', 'Residual'),2),
                    out_of = factor(c(rep('total', 4), rep('RE', 4)), levels = c('total', 'RE')),
                    percent = c(c(id_v, idsess_v, other_RE, s_2)*100 / total,
                                c(id_v, idsess_v, other_RE, NA)*100 / total_RE))
  if(other_RE == 0){
    rez <- rez[rez$source != 'Other RE',]
  }
  return(rez)
})

The plots on the left, below, show model-expected functional connectivity pearson correlations for each participant, for each session, overlayed on the raw data (one point per parcel, per session, per participant). A line for each participant's model-expected mean across all sessions is overlayed on this. The right plot shows proportions of variance accounted for by the random effect estimates as a proportion of both total variance, and total random effects (RE) variance. The total RE variance includes individual means (ID), individual means per session (ID/Session), and the difference in connectivity between right and left hemispheres (we treat this as a nuisance variable, essentially).

library(patchwork)

plot_percents <- function(percent_table){
  ggplot(percent_table, aes(y = percent, x = out_of, fill = source)) + 
    geom_col(position = 'stack') + 
    scale_fill_manual(breaks = c('ID', 'ID/Session', 'Other RE', 'Residual'),
                      values = apal[c(4,3,2,1)]) +
    scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100), labels = paste0(c(0, 20, 40, 60, 80, 100), '%')) + 
    labs(x = 'Variance (denominator)', y = '', fill = 'Variance source') +
    jftheme
}

plot_predictions <- function(amodel, point_alpha){
  adt <- as.data.table(amodel@frame)
  adt[, wave := factor(as.numeric(gsub('ses-(\\d+)', '\\1', sess)))]
  newdata <- unique(adt, by = c('id', 'sess', 'wave'))[, c('H', 'z') := list(0, NULL)]
  newdata$z_prime <- predict(amodel, newdata = newdata)
  newdata$z_prime_id <- predict(amodel, newdata = newdata, re.form = ~(1|id))
  ggplot(newdata, aes(x = wave, y = tanh(z_prime), group = id)) + 
    #geom_violin(data = adt, aes(group = NULL, y = tanh(z)), size = 0, alpha = .5, fill = apal[[1]]) +
    geom_point(data = adt, aes(group = NULL, y = tanh(z)), alpha = point_alpha, color = apal[[1]], position = position_jitter(), size = .5) +
    geom_point(color = apal[[3]], alpha = 1) + 
    geom_line(color = apal[[3]]) + 
    geom_rug(data = unique(newdata, by = 'id'), aes(y = z_prime_id, x = NULL), color = apal[[4]], alpha = 1, sides = 'tr', length = unit(0.03, "npc")) + 
    geom_hline(yintercept = 0, color = apal[[1]]) + 
    coord_cartesian(ylim = c(-.025, .65), xlim = c(1, 10)) + 
    labs(y = 'FC correlation', x = 'Session') + 
    jftheme
}

max_points <- unlist(sea_schaefer400_7_withinnet[, .N, by = net1][, list(max = max(N))])
point_alphas <- sea_schaefer400_7_withinnet[, .N, by = net1][, p := 1 - N / max_points][, pomp := (p - min(p)) / (max(p) - min(p))][, alpha := .025 + pomp*.1]

font_add_google("Didact Gothic", "Didact Gothic")
showtext_auto()
for(i in 1:length(model_fits)){
  model_fit <- model_fits[[i]]
  network_name <- networks[[i]]
  cat('\n\n### ', network_name, '{.tabset}\n\n') 
  cat('\n\n#### Plot\n\n')
  point_alpha <- point_alphas[net1 == network_name, alpha]
  print(plot_predictions(model_fit, point_alpha = point_alpha) + plot_percents(proportion_variance_tables[[i]]) +
          plot_layout(ncol = 2, widths = c(4,1)))
  cat('\n\n#### Table (%)\n\n')
  print(knitr::kable(tidyr::spread(proportion_variance_tables[[i]], out_of, percent), digits = 1))
  cat('\n\n#### Model Summary\n\n')
  cat('\n```\n')
  print(summary(model_fit))
  cat('\n```\n')
}

Cont

Plot

Table (%)

source	total	RE
ID	0.0	0
ID/Session	6.6	100
Residual	93.4	NA

Model Summary

Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
   Data: this_net_dt

REML criterion at convergence: 48899.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.5451 -0.6778 -0.0803  0.5986  5.6264 

Random effects:
 Groups   Name        Variance Std.Dev.
 sess:id  (Intercept) 0.005433 0.07371 
 id       (Intercept) 0.000000 0.00000 
 Residual             0.076803 0.27713 
Number of obs: 176231, groups:  sess:id, 281; id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.196543   0.004447   44.19
convergence code: 0
boundary (singular) fit: see ?isSingular

Default

Plot

Table (%)

source	total	RE
ID	0.0	0
ID/Session	5.9	100
Residual	94.1	NA

Model Summary

Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
   Data: this_net_dt

REML criterion at convergence: 181397.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.7510 -0.6822 -0.0861  0.5938  6.3337 

Random effects:
 Groups   Name        Variance Std.Dev.
 sess:id  (Intercept) 0.005088 0.07133 
 id       (Intercept) 0.000000 0.00000 
 Residual             0.081313 0.28516 
Number of obs: 548180, groups:  sess:id, 281; id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.185122   0.004273   43.32
convergence code: 0
boundary (singular) fit: see ?isSingular

DorsAttn

Plot

Table (%)

source	total	RE
ID	0.1	1.3
ID/Session	8.7	98.7
Residual	91.2	NA

Model Summary

Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
   Data: this_net_dt

REML criterion at convergence: 60794.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.9172 -0.6943 -0.0936  0.5954  5.6176 

Random effects:
 Groups   Name        Variance  Std.Dev.
 sess:id  (Intercept) 0.0085855 0.09266 
 id       (Intercept) 0.0001118 0.01058 
 Residual             0.0904365 0.30073 
Number of obs: 137318, groups:  sess:id, 281; id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.250908   0.005914   42.42

Limbic

Plot

Table (%)

source	total	RE
ID	0.0	0
ID/Session	13.1	100
Residual	86.9	NA

Model Summary

Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
   Data: this_net_dt

REML criterion at convergence: -1918.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0608 -0.6610 -0.0535  0.5952  5.3055 

Random effects:
 Groups   Name        Variance Std.Dev.
 sess:id  (Intercept) 0.008228 0.09071 
 id       (Intercept) 0.000000 0.00000 
 Residual             0.054555 0.23357 
Number of obs: 38920, groups:  sess:id, 266; id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.189825   0.005714   33.22
convergence code: 0
boundary (singular) fit: see ?isSingular

SalVentAttn

Plot

Table (%)

source	total	RE
ID	0.0	0
ID/Session	8.8	100
Residual	91.2	NA

Model Summary

Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
   Data: this_net_dt

REML criterion at convergence: 27433.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.9508 -0.6711 -0.0699  0.6001  6.6801 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 sess:id  (Intercept) 6.739e-03 8.209e-02
 id       (Intercept) 7.511e-11 8.667e-06
 Residual             7.019e-02 2.649e-01
Number of obs: 145202, groups:  sess:id, 281; id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.221196   0.004947   44.72
convergence code: 0
boundary (singular) fit: see ?isSingular

SomMot

Plot

Table (%)

source	total	RE
ID	0.0	0
ID/Session	10.6	100
Residual	89.4	NA

Model Summary

Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
   Data: this_net_dt

REML criterion at convergence: 132084.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.2412 -0.6741 -0.1112  0.5556  7.7342 

Random effects:
 Groups   Name        Variance  Std.Dev. 
 sess:id  (Intercept) 9.641e-03 9.819e-02
 id       (Intercept) 8.132e-11 9.018e-06
 Residual             8.170e-02 2.858e-01
Number of obs: 392102, groups:  sess:id, 281; id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.278110   0.005875   47.34
convergence code: 0
boundary (singular) fit: see ?isSingular

Vis

Plot

Table (%)

source	total	RE
ID	0.0	0
ID/Session	4.7	100
Residual	95.3	NA

Model Summary

Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
   Data: this_net_dt

REML criterion at convergence: 169074.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7379 -0.7070 -0.1085  0.6245  5.5535 

Random effects:
 Groups   Name        Variance Std.Dev.
 sess:id  (Intercept) 0.005786 0.07606 
 id       (Intercept) 0.000000 0.00000 
 Residual             0.117966 0.34346 
Number of obs: 239843, groups:  sess:id, 281; id, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.268284   0.004594    58.4
convergence code: 0
boundary (singular) fit: see ?isSingular